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We propose a chemical potential dependent effective gluon propagator and study 
the chiral quark condensate in strongly interacting matter in the framework of 
Dyson-Schwinger equation approach. The obtained results manifest that, as the 
effect of the chemical potential on the effective gluon propagator is taken into ac- 
count, the chiral quark condensate decreases gradually with the increasing of the 
chemical potential if it is less than the critical value, and the condensate vanishes 
suddenly at the critical chemical potential. The inclusion of the chemical potential 
in the effective gluon propagator enhances the decreasing rate and decreases the 
critical chemical potential. It indicates that the chiral symmetry can be restored 
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completely at a critical chemical potential and restored partially as the chemical 
potential is less than the critical value. If the effective gluon propagator is indepen- 
dent of the chemical potential, the chiral symmetry can only be restored suddenly 
but no gradual restoration. 

PACS numbers: 24.85. +p, 12.38.Aw, 12.38.Lg, 03.75.Hh 



2 



The chiral symmetry spontaneous breaking and color confinement have been known as 
two essential characters of strong interaction in low energy region. The relation between 
these two essential properties is a fundamental problem in strong interaction physics. 
However, it is still very difficult to study such a relation directly at present. On the 
other hand, it is fortunate that one has known the vacuum of strong interaction would 
become trivial at high enough temperature and/or chemical potential because of the 
very weak-coupling interaction, or in other word, the asymptotic freedom. Then the 
behaviors of the chiral symmetry restoration and color deconfinement at hight temperature 
and/or chemical potential would give signatures or information about the relation between 
the two properties. Concerning the measurement of the chiral symmetry restoration, 
one used to take the chiral quark condensate (qq) as a characteristic since it is believed 
to be one of the most important configuration of the strong interaction vacuum (see 
for example Refs.[H [2]). As for the methodology, the Lattice QCD simulation can be 
safely expand to finite temperature [31 H] and has made great efforts to deal with the 
problem at finite chemical potential[5l El El O EH [11]. With the difference between the 
fundamental and adjoint representations of the color symmetry of a quark being taken into 
account, lattice QCD simulation and SU(3) gauge theory analysis show that the quark 
deconfinement and chiral symmetry restoration can happen at same critical temperature 
or not[3l[T2]. Meanwhile many approaches in the framework of continuous field theory of 
QCD have been implemented to study the chiral symmetry breaking and its restoration in 
dense strongly interacting matter at zero temperature, and different models gave different 
behaviors at finite chemical potential (see for example Refs. (13], [HI [151 EH EL71 [HI [HI 1201 
1211 [221 l23"l [2H [251 [26]). However, there is still not any work to take into account the effect 
of the chemical potential dependence of the effective gluon propagator. 

It has been shown that the Dyson-Schwinger (D-S) equation approach provides a non- 
perturbative framework which admits a simultaneous study of dynamical chiral symmetry 
breaking and color confinement |27[ [28] . With nonperturbative truncations preserving the 
chiral symmetry, this approach has been widely used to study the properties of strong 
interaction vacuum and the properties of hadrons in free space [29] [30] . Meanwhile it has 
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been extended to the system with finite temperature and/or finite chemical potential to 
simulate the chiral symmetry restoration and deconfinement [311 132] . In this paper we will 
then take the Dyson- Schwinger equation approach at finite chemical potential to calculate 
the chemical potential dependence of the chiral quark condensate in strongly interacting 
matter. In previous studies on the chiral behavior at finite chemical potential, the ef- 
fective gluon propagator is usually taken as the same as that in free space. We should 
know that the chemical potential does infect the interaction. Then the chemical potential 
dependence of the gluon propagator has to be deliberated. In this paper we propose an 
approximation to take into account the chemical potential dependence of the effective 
gluon propagator and investigate the effect of the interaction strength and the screening 
parameters in the modified effective gluon propagator on the chiral quark condensate at 
finite chemical potential. 

As the lowest dimension condensate of quarks and gluons, the chiral quark condensate 
is the essential characteristic in the QCD phase transition, and hence plays essential role 
in describing hadron structure and properties of nuclear matter and finite nuclei. In the 
chiral limit, the chiral quark condensate can be evaluated from the quark mass function 
in ultraviolet region at zero chemical potential or with the definition 



where G(q) is the dressed quark propagator at zero chemical potential. It has been 
shown that, with the renormalized D-S equation, such an expression is equivalent to 
the one deduced from the mass function in QCD Sum rules [33]. In the case of finite 
chemical potential, one usually takes the same expression to study the chemical potential 
dependence of the chiral quark condensate [20], 12H [33]. We have then 



where Q[fj](q) is the dressed quark propagator at finite chemical potential. It is then 
imperative to investigate the chemical potential dependence of the dressed quark propa- 
gator. 




(1) 
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The main aim of D-S equation approach is just to study the propagators of fermion and 
boson at all energy scale. The D-S equation approach is based on a coupled set of integral 
equations among quark, gluon, ghost and vertex functions. They form a countably infinite 
set of coupled integral equations with the one for a n-point Schwinger function depending 
on the (n + 1) and higher point functions. In order to handle the system practically, it is 
necessary to make certain simplifications and truncations. One simple approach, which 
is commonly referred to as the rainbow approximation, employs the bare quark-gluon 
vertex and leaves the equation for the dressed quark propagator to be solved with a given 
effective gluon propagator as input. Under the rainbow approximation, the dressed quark 
propagator at zero chemical potential G(p) = = 0](p) can be determined well with 
the truncated D-S equation [511 ESI EE] 



where D^k) is an effective gluon propagator, and the inverse of the dressed quark prop- 
agator is conventionally decomposed as 



with A(p 2 ) and B(p 2 ) being scalar functions. The chiral quark condensate at zero chemical 
potential can be determined by the two scalar functions as: 



Basing on the solution of the D-S equation, one can obtain a trivial solution B(p 2 ) = 0, 
namely the Wigner solution that characterizes a phase in which chiral symmetry is not 
broken and the dressed quarks are not confined. This phase relates to the trivial vacuum 
with (qq) = 0. We are mainly interested in the nontrivial solution with nonzero 5(Nambu 
solution) which should be found by numerical iterations. This solution corresponds to 
the nontrivial vacuum containing the dynamical quark mass generation, the chiral quark 
condensate and the existence of massless Goldstone bosons. It is also necessary to mention 
that the Wigner solution is always possible due to the path-integral formula used in the 
framework of the D-S equation. Nevertheless, the Nambu solution is only possible if 






(4) 
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the coupling is strong enough at the infrared region [271 (37J [38]. Since the chiral quark 
condensate is zero in Wigner phase and nonzero in Nambu phase, the variation behavior 
of the chiral quark condensate with respect to the chemical potential can simulate the 
chiral phase transition from the Nambu phase to Wigner phase. 

The chemical potential is introduced as a Lagrange multiplier pxpi^q in QCD action. 
Similarly, the dressed quark propagator G\p](p) at non-zero chemical potential fi can be 
written as 



4 f d 4 q 

3~ MO) = *7 • v - 74/i + - J j^z D ^(p - QhnSbAteh* 



(6) 



where the form of the effective interaction D^ v is usually taken as the same as that in 
Eq. ©. 

In prevenient approach of D-S equation at finite chemical potential, the dressed-quark 
propagator can be written, in general, as[20j, [31] 

Q~ X Wp) = %l ■ pA(p) + i lA (p A + ifi)C(p) + B(p) , (7) 

where p u = (p,p A + i/x), the complex functions A, C and B can be defined by the quark 
equation with an effective gluon propagator. With a 5-function in momentum space for 
the effective gluon propagator (Munczek-Nemirovsky model [2H]), h has been found that, 
in the D-S equation approach [201 [24] . the quark condensate increases with chemical 
potential, which is not consistent with the behavior given in composite operator approach 
of QCD [13J, dilute instanton liquid model [2JJ, Nanbu-Jona-Lasinio (NJL) model [23] 
and other approaches [TH [TBI E21 [25]. Looking over the characteristic of the Munczek- 
Nomirovsky model, which gives an infinity at zero exchanged momentum, we propose 
that the discrepancy of the variation behaviors of the chiral quark condensate against the 
chemical potential arises from the overestimated sudden enhancement in the far infrared 
region. To remove the discrepancy, we should then suppress the far infrared enhancement. 
Prevenient works have shown that the model with effective gluon propagator 



D^(k)=t,Av 2 d-^-^, (8) 

with d = and t^ v = 5^, which means that the "Feynman-like" gauge is taken in prac- 
tical calculation, can describe the pion weak decay constant and other low energy chiral 



observables well (see for example Ref.jlO]). It is obvious that this effective gluon propa- 
gator has an finite infrared enhancement and does not involve the ultraviolet behavior of 
the QCD running coupling. Since such a model produces ultraviolet convergent integrals 
naturally, the renormalization is not necessary. Meanwhile, this effective gluon propaga- 
tor does not have a Lehmann representation, and the corresponding classical potential 
between quarks can be written as 



It is apparent that such a potential corresponds to an approximately quadratic con- 
finement in intermediate range (up to r « 3.0 fm) and linear-like in very short range 
(r < 0.2 fm). In most recent years, great progress has been made in systematic studies 
on quark-quark interaction (see for example Ref. [3U [351 HD H21 H31 HU H5] ) . One should, 
in principle, take the obtained analytic properties of the gluon and quark propagators 
(c.f. given in Ref. [3^1 135]) to study the chiral phase transition. However, for simplicity 
in studying the quark D-S equation and incorporating the chemical potential effect at 
present stage, we take the effective gluon propagator in Eq. (jHJ) as the starting point of 
our model of interaction. 

In general approach to study the quark equation at finite chemical potential, one 
should take a chemical potential \i dependent effective gluon propagator as input, since 
the chemical potential of the matter also influences the interaction. However, the related 
knowledge is still in lack up to now. For simplicity, we propose an approximation for the 
effective gluon propagator at finite chemical potential with an extension k 2 — > k 2 + /3/j, 2 , 
the Eq. (JBJ) is then modified as 



where /i denotes the chemical potential and (3 is a scaling parameter denoting the strength 
of the /i dependence. The classical potential with this modified gluon propagator can be 
written as 




(9) 



D'{k) = t^ 2 d 




(10) 





7 



It is evident that the strength of the attractive interaction becomes smaller with an 
increase of the chemical potential and the scaling parameter (3. 

Substituting Eq. ffTUl) into Eq. and accomplishing some derivations, we have the 
equations for the scalar functions A(p), B{p) and C{p) in Eq. ([7]) as 

(A(v) - l)r? = - I d " q ^ 2xH Am ' t (12) 

1 {P) )P 3 J (2 7 r) 4 (A; 2 + / 9/i 2 ) 2 + A ( fA 2 (g) + (g4 + ^) 2 C 2 (g) + J B 2 (g)' 1 } 



(r(~\-\M ■ \ - 8 / ^ ^ 2 X 2 d C(q)(q A + ifi) 

[ [P) L)[P4 + W 3 J (27r) 4 (fc 2 +^ 2 ) 2 + Ag^(g) + (g 4 +^)2C 2 (g) + J B 2 (g)' 1 ' 

yP> 3 J (27r) 4 (P + /3/i 2 ) 2 + Ag 2 A 2 (g) + (g 4 + z/i) 2 C 2 (g) + J B 2 (g) 1 } 
With the parameter being taken as A = 0.01 GeV 4 , x = 1-33 GeV, with which the pion 
decay constant can be fitted to 87 MeV (quite close to the experimental date 93 MeV) [4"0] . 
we can solve the coupled equations at zero chemical potential easily and obtain the chiral 
quark condensate in vacuum as (250 MeV) 3 . It is evident that such a result agrees 
excellently well with the empirical value. 

At nonzero chemical potential, it is a quite hard task to solve the coupled integral 
equations (fT2"|) - fTH|) . since they depend on the momenta p, the chemical potential fi as 



well as the angle 6 between p and fi with cos^ = p ■ ^/ \/p 2 jJ? due to the breaking of 
the 0(4) symmetry in the four momentum space. By abandoning the commonly used 
discretion of the unknown functions but implementing the smooth polynomial approxima- 
tions^, 46J for the to be determined complex functions A(p), C(p) and B(p), we solve 
the equations f|T2l - f|T4l) (for the process of the smoothing and solving the equation, please 
see the Appendix-A). 

With the solutions of the D-S equations at any chemical potential and scaling param- 
eter p, we evaluate the chemical potential dependence of the chiral quark condensate at a 
series values of the scaling parameter (3. A part of the obtained results are illustrated in 
Fig. [TJ The figure shows apparently that, as the chemical potential increases in the region 
less than a critical value, the in-medium chiral quark condensate decreases monotonously 
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Figure 1: Calculated chemical potential dependence of the chiral quark condensate at 
several scaling parameters 

for any value of (3 6 (0, 1]. As the critical value of the chemical potential /i c is reached, 
the chiral quark condensate vanishes suddenly. Such a behavior indicates that the chiral 
symmetry can be restored partially before the critical chemical potential is reached. At 
the critical point of the chemical potential, the chiral symmetry can be completely re- 
stored suddenly. And the phase transition is in first order. Such a result is consistent with 
that given in the composite operator formalism of QCD [13] . Moreover, with the increase 
of the scaling parameter (3, the decreasing rate of the chiral quark condensate against 
the chemical potential gets obviously larger and the critical chemical potential becomes 
apparently smaller. For instance, for (3 = 0.3, 1.0, the critical chemical potential /i c is 
344 MeV, 245 MeV and the ratio | (qq)^\/\(qq)o\ is 72.4%, 56.3%, respectively. In addition, 
in the case of that the chemical potential does not affect the effective gluon propagator, 
i.e., (3 = 0, the chiral quark condensate in the matter maintains almost a constant if the 
chemical potential is less than the critical value (/i c = 546 MeV, \(qq) /j, c \/\(qq)o\ — 99.6%). 
Such a result is analogous to that given in the NJL model [23] and dilute instanton liquid 
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model [21] - These results indicate that the critical chemical potential for the chiral quark 
condensate to be restored completely decreases with the increase of the scaling parame- 
ter f3. The obtained result of the scaling parameter dependence of the critical chemical 
potential is illustrated in Fig. [2J 




Figure 2: The scaling parameter dependence of the critical chemical potential for the 
chiral symmetry to be restored 

To understand the variation characteristics of the chiral quark condensate against the 
chemical potential and the scaling parameter intuitively, we recall the definition of the 
condensate in Eq. (j2J). It is evident that the integrand for the condensate reads 

^A^{p) + { Pi + i^yC\p)+B^{p)' K ] 

where the A(p), C(p) and B{p) are the solutions of the D-S equations (IT2l - (IT%j) . After 
integrating over the angle, one can express the momentum dependence of the integrand 
explicitly as (neglecting some constants) 

IG(f) = 1*0 sin* ^ A2m + (p f + B i ^ c2(]}) + BHij) , (16) 
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with cos 6 = p- fi/ ^/p 2 fi 2 . Since the real part of both A(p) and C(p) is even and increasing 
function of the cos#, the real part of B(p) is an even but decreasing function of the cos#, 
and the imaginary part of A(p), C(p) and B(p) is odd function of cos#, such an integrand 
is a real function of the momentum p. The behavior of the function IG(p 2 ) at several 
values of the chemical potential and the scaling parameter can be displayed in Fig. [31 It is 




Figure 3: Calculated chemical potential dependence of IG(p 2 ) at several scaling parame- 
ters: solid curve for (3 = 0.0; dash curve for (3 = 0.3; dot curve for (3 = 1.0. 

evident that, at any chemical potential, the function IG(p 2 ) is convergent in the infrared 
and the ultraviolet regions. It indicates that the condensate is well defined. In more detail, 
there exists a critical chemical potential /z c , if /i > fi c , IG{p 2 ) = 0. Such a fi c decreases 
with the increasing of the scaling parameter (3, for example, /i c = 0.546 GeV for (3 = 0.0, 
// c = 0.344 GeV for (3 = 0.3, /i c = 0.245 GeV for (3 = 1.0. In the case of (3 + 0, as /i < /i c , 
IG(p 2 ) increases with the increasing of momentum. However the IG(p 2 ) decreases at large 
momentum. The momentum for the IG(p 2 ) to take the maximal value shifts to smaller 
as the chemical potential and the scaling parameter increase. Moreover, the maximal 
value of IG(p 2 ) decreases as the scaling parameter increases. These characteristics of 
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the integrand manifests evidently that, if the chemical potential is smaller than a critical 
value but influences the effective gluon propagator (i.e., with (3 G (0, 1]), the chiral quark 
condensate decreases with the increasing of chemical potential, and the decreasing rate 
gets larger as the scaling parameter increases. As the chemical potential reaches the 
critical value, the chiral quark condensate vanishes suddenly. In the case of that the 
chemical potential does not influence the effective gluon propagator (i.e., with (3 = 0), 
even though the solutions of the D-S equation at finite chemical potential are different 
from those at zero chemical potential, because of the compensation of the solutions A, C 
and B at the presently chosen parameters, the integrand of the chiral quark condensate 
behaves almost independent of the chemical potential if it is less than the critical value. 
The chiral quark condensate turns to be then approximately a constant before it vanishes 
suddenly at the critical chemical potential. On the other hand, concerning the mass 
function of a quark in the case of that the chemical potential does not affect the effective 
gluon propagator, we realize that it does not change drastically either with the chemical 
potential as it is smaller than the critical value. It means that the dressing effect come 
from the chiral quark condensate does not change obviously. It hints then the chiral quark 
condensate maintains almost a constant if the effective gluon propagator is independent 
of the chemical potential. 

To explore the underlying physics of the variation behavior of the chiral quark con- 
densate against the chemical potential and the scaling parameter, we look through the 
characteristic of the model effective gluon propagator. From Eqs. ( 1101) and fill I) , one can 
recognize easily that the model effective gluon propagator is the one with approximately 
quadratic confinement in intermediate region and linear-like in very short region. Both 
the chemical potential and the scaling parameter play roles of screening on the confine- 
ment. It means that the increases of the chemical potential and the scaling parameter 
decrease the confinement strength. It has been known that the chiral quark condensate 
is proportional to the mean spectral density at zero energy of quarks with Banks-Casher 
relation [33, HZ] ~(qq) — 7ri/ g (0) . In general principle, the mean spectral density i/(0) de- 
creases as the attractive interaction gets weaker. The classical potential in Eq. ffTTj) shows 
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that, in the case of that the effective gluon propagator depends on the chemical potential 
(i.e., with scaling parameter (3 > 0), the interaction in the system is just attractive and 
gets weaker as the chemical potential and the scaling parameter increase. Corresponding 
to the increases of the chemical potential and the scaling parameter, the mean spectral 
density at zero energy of quarks decreases, so that the absolute value of the chiral quark 
condensate decreases. Moreover, as the chemical potential is large enough (i.e., the den- 
sity of quarks is sufficiently large), the screening effect is so strong that the interaction 
between quarks becomes very weak, even vanishes (asymptotic free), which is analogous 
to the molecular in a liquid, where it experiences an attraction to the inner if it is in the 
surface region and becomes almost free in the inner region. As a consequence, the chiral 
quark condensate changes to zero suddenly. Combining the calculated variation behavior 
of the chiral quark condensate with respect to the chemical potential (or the density of 
the strong interaction matter) and the inferred changing characteristic of the interaction 
between quarks, we can recognize that the chiral symmetry can be restored abruptly as 
the density of the matter reaches a critical value, and the chiral phase transition is in first 
order. In addition, the chiral symmetry can be restored partially and gradually before the 
phase transition takes place. However, in the case of that the effective gluon propagator 
is independent of the chemical potential (i.e., with /3 = 0), the interaction is free from the 
chemical potential, or in other word, the chemical potential does not screen the attractive 
interaction. Then the chiral quark condensate almost maintains a constant before the 
critical chemical potential is reached. 

In summary, we have calculated the chiral quark condensates in vacuum and in 
strongly interacting matter in the Dyson-Schwinger equation approach with a model 
effective gluon propagator including finite chemical potential effect in this paper. The 
calculated results show that, as the chemical potential influences the effective gluon prop- 
agator, the chiral quark condensate decreases with the increase of the chemical potential 
before the critical value is reached and the condensate vanishes suddenly at the critical 
chemical potential. Moreover, the increase of the scaling parameter of the chemical po- 
tential in the effective gluon propagator enhances the decreasing rate and decreases the 
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critical chemical potential. It indicates that the chiral symmetry can be restored abruptly 
as the density of the matter reaches a critical value, and the chiral phase transition is in 
first order. Before the critical density or chemical potential is reached, the chiral symme- 
try can be restored partially and gradually. In addition, as the chemical potential does not 
affect the effective gluon propagator, the chiral symmetry can only be restored suddenly 
at a critical chemical potential, but no gradual restoration takes place before the critical 
chemical potential is reached. It manifests that the effect of the chemical potential on 
the gluon propagator plays significant role in the process of the chiral symmetry restora- 
tion. Concerning the modification on the effective gluon propagator at finite chemical 
potential, we just extend the exchanged momentum related term from k 2 to chemical 
potential dependent k 2 + j3fi 2 . The parameter f3 or the Pfj, 2 displays the effect of the finite 
chemical potential on the effective gluon propagator or the quark-quark interaction. This 
is definitely only a model or an approximation. Furthermore, we have taken the rainbow 
approximation for the quark equation. Then studying the strong interaction matter by 
solving the coupled quark, gluon and ghost equations systematically is imperative. The 
related investigations are under progress. 

This work was supported by the National Natural Science Foundation of China 
(NSFC) under contract Nos. 10425521 and 10575004, the Major State Basic Research 
Development Program under contract No. G2000077400, the research foundation of the 
Ministry of Education, China (MOEC), under contact No. 305001 and the Research Fund 
for the Doctoral Program of Higher Education of China under grant No 20040001010. 
One of the authors (YXL) thanks also the support of the Foundation for University Key 
Teacher by the MOEC. 

Appendix-A: 

Because the 0(4) symmetry in the four momentum is broken, the complex functions 
A(p), C(p) and B{p) depend not only on the square of momentum p 2 and chemical 
potential /x 2 but also on the angle 6 between p and jj with cos^ = p ■ p,j ^Jp 2 ^ 2 . To solve 
the coupled equations (|T2|) - (1T4|) . a smooth polynomial approximation should be taken for 
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both the angle and the square of momentum when we discrete the integrand in numerical 
calculation. 

The angle 9 dependence of the functions is analyzed by an expansion in terms of the 
second kind Chebyshev polynomials C/j(cos0) as 

N ch>2 -i 

A(p 2 )= A m (p 2 ;V 2 )U m (cos6)(i) m , (Al) 



m=0 
N ch , 2 -1 



C(P 2 )= E C m {p 2 -^ 2 )U m {cos6)(iT, (A2) 

N ch , 2 -1 

B(p 2 ) = E B rn {p 2 ^ 2 )U m {cos6)(i) m . (A3) 

Since the Uk(cos8) form an orthogonal set, we can project the equations (|T2|) . (fT3|) and 
(THI) onto a set of nonlinear integral equations satisfied by the form factors A m (p 2 ; fi 2 ), 
C m (p 2 ; /i 2 ) and B m (p 2 ; fi 2 ) which are all real functions depending only on p 2 and /i 2 . In 
practical calculation we truncate the expansion at a certain rank N c h^ and the degree of 
the fitting can be mitigated by increasing N^^- After the angular expansion we then take 
Chebyshev expansion for the form factors A m (p 2 ; fi 2 ), C m (p 2 ;fi 2 ) and B m (p 2 ; fi 2 ), which 
gives 



Ami?'; = E <TMP)) - ( A4 ) 

N chA x 

C m {p 2 ^ 2 ) = E - Y' ( A5 ) 

j'=i 

with 

log(p 2 /Ae) 



s(p) 



log(A/e) ' 

where Tj(s(p)) is the j-rank Chebyshev polynomial of the first kind, A is the ultraviolet 
cutoff, and e is the infrared cutoff. In this expansion we make use of the first kind 
Chebyshev polynomial with a cutoff N^i. Up to now the unknown functions A, C and 
B have been expressed by the smooth polynomial with 3N c h y iN c h,2 Chebyshev coefficients 
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a J m , c J m and b> m which would be numerically calculated by an effective iteration method. 
In practical calculation, we take N^i = 32, = 6, A = 10 GeV and e = 0.01 GeV 

(considering the parameter A = O.OlGeV 4 and e 4 = 10~ 8 GeV 4 , such a choice of the cutoff 
is quite reasonable). 

In order to determine the SN c h t iN c h t 2 Chebyshev coefficients a 3 m , c J m and b 3 m , we require 
these integral equations to satisfy N ch; i fixed external momenta. To solve the nonlinear 
equations, we implement the Newton method. This method takes derivatives of the equa- 
tions with respect to the unknowns to speed up the convergence. After some derivations, 
it is transformed to a set of ((3N ch:1 N ch:2 ) x {3N ch ^N ch)2 )) linear equations, which can be 
solved easily. As for the integration we introduce spherical coordinates, i.e. the volume 
element d 4 q reads 2nq 3 dq sin 2 9d9 sin <pd(f). The integration ranges of the variables are: 
q G [e, A] and 9, <fi G [0,7r]. The integration is calculated by the N-points Gauss- Legendre 
quadrature rule. The integrand may still not be smooth enough at the junctions of some 
intervals, it is then necessary to split the region for much higher accuracy. We split the in- 
tegrand into three regions [e 2 , min(p 2 , q 2 )}, [min(p 2 , q 2 ),max(p 2 , q 2 )} and [max(p 2 , q 2 ), A 2 ], 
where p 2 is the square of the external momentum and q 2 is defined by the possible zero 
of q*A 2 (q) + (q 4 + ifi) 2 C 2 (q) + B 2 (q) at q 4 = 0. 
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